/*
 *  MWRK.cpp
 *  T3nsors2
 *
 *  Created by Michael Barriault on 10-06-30.
 *  Copyright 2010 University of Guelph. All rights reserved.
 *
 */

#include "MW.h"
#include <iostream>
using namespace t3;
using std::cout;
using std::endl;

void MW::RK(double dt) {
	scalar *j1 = PhiDot(Phi, Psi, E, B);
	scalar *k1 = PsiDot(Phi, Psi, E, B);
	vector *l1 = EDot(Phi, Psi, E, B);
	vector *m1 = BDot(Phi, Psi, E, B);

	Bounds(j1,k1,l1,m1,true,true,false);

	scalar *Phi2 = scalar::init(O,"Phi2");
	scalar *Psi2 = scalar::init(O,"Psi2");
	vector *E2 = vector::init(O,3,"E2");
	vector *B2 = vector::init(O,3,"B2");
	FOR(o,O->size()) Phi2->at(o) = Phi->at(o) + 0.5*dt*j1->at(o);
	FOR(o,O->size()) Psi2->at(o) = Psi->at(o) + 0.5*dt*k1->at(o);
	FOR(a,3) FOR(o,O->size()) E2->c(a)->at(o) = E->c(a)->at(o) + 0.5*dt*l1->c(a)->at(o);
	FOR(a,3) FOR(o,O->size()) B2->c(a)->at(o) = B->c(a)->at(o) + 0.5*dt*m1->c(a)->at(o);

	scalar* j2 = PhiDot(Phi2, Psi2, E2, B2);
	scalar* k2 = PsiDot(Phi2, Psi2, E2, B2);
	vector* l2 = EDot(Phi2, Psi2, E2, B2);
	vector* m2 = BDot(Phi2, Psi2, E2, B2);
	
	Bounds(j2,k2,l2,m2,true,true,false);

	scalar *Phi3 = scalar::init(O,"Phi3");
	scalar *Psi3 = scalar::init(O,"Psi3");
	vector *E3 = vector::init(O,3,"E3");
	vector *B3 = vector::init(O,3,"B3");
	FOR(o,O->size()) Phi3->at(o) = Phi->at(o) + 0.75*dt*j2->at(o);
	FOR(o,O->size()) Psi3->at(o) = Psi->at(o) + 0.75*dt*k2->at(o);
	FOR(a,3) FOR(o,O->size()) E3->c(a)->at(o) = E->c(a)->at(o) + 0.75*dt*l2->c(a)->at(o);
	FOR(a,3) FOR(o,O->size()) B3->c(a)->at(o) = B->c(a)->at(o) + 0.75*dt*m2->c(a)->at(o);
	
	scalar* j3 = PhiDot(Phi3, Psi3, E3, B3);
	scalar* k3 = PsiDot(Phi3, Psi3, E3, B3);
	vector* l3 = EDot(Phi3, Psi3, E3, B3);
	vector* m3 = BDot(Phi3, Psi3, E3, B3);
	
	Bounds(j3,k3,l3,m3,true,true,false);

	FOR(o,O->size()) Phi->at(o) += dt*(2*j1->at(o)+3*j2->at(o)+4*j3->at(o))/9;
	FOR(o,O->size()) Psi->at(o) += dt*(2*k1->at(o)+3*k2->at(o)+4*k3->at(o))/9;
	FOR(a,3) FOR(o,O->size()) E->c(a)->at(o) += dt*(2*l1->c(a)->at(o)+3*l2->c(a)->at(o)+4*l3->c(a)->at(o))/9;
	FOR(a,3) FOR(o,O->size()) B->c(a)->at(o) += dt*(2*m1->c(a)->at(o)+3*m2->c(a)->at(o)+4*m3->c(a)->at(o))/9;
	
	j1->release(); j2->release(); j3->release();
	k1->release(); k2->release(); k3->release();
	l1->release(); l2->release(); l3->release();
	m1->release(); m2->release(); m3->release();
	Phi2->release(); Phi3->release();
	Psi2->release(); Psi3->release();
	E2->release(); E3->release();
	B2->release(); B3->release();
}
